; test which roots are the physical solutions 
;==========================================================


; regular (non tau0) analytic solution for xHI
; it seems that there is always one solution between 
; 0 and 1 and another above 1 (can I prove this)
;==========================================================
GHI = 1.0d-12
T = 1.0d2
y = 0.0d0
nH = 1.0d-1

RC = hiirca_hui(T)
CI = hici_hui(T)

R = (CI + RC) * nH
Q = -( GHI + (CI+2*RC) * nH + (CI+RC) * nH * y )
P = RC * nH * (1.0d0 + y)
d = Q*Q - 4 * P * R

xHIp = (-Q + sqrt(d) ) / (2*R)
xHIm = (-Q - sqrt(d) ) / (2*R)   ; NEGATIVE ROOT! 

print
print
print, '-----------------------------------------'
print, 'case of regular solution'
print
print, 'log xHIp = ', alog10(xHIp)
print, 'log xHIm = ', alog10(xHIm)
print

qq = - 0.5d0 * (Q + abs(Q)/Q * sqrt(d))

xHIp = qq / R
xHIm = P / qq   ; NEGATIVE ROOT! 

print, 'log xHIp = ', alog10(xHIp)
print, 'log xHIm = ', alog10(xHIm)
print
print, '-----------------------------------------'


; plot zeros of non tau0 solution for several photoionization rates
;----------------------------------------------------------
i_plot_multi=1

GHI_arr = [-23.0d0, -21.0d0,-14.0d0, -13.6d0, -13.4d0, -13.2d0]
GHI_arr = [-22.0d0, -21.0d0, -14.0d0, -13.0d0,-12.0d0,-11.0d0]
GHI_arr = [-12.0d0]

log_T_min = 2.0d0
d_log_T = 0.1d0

N=1000000
Tmodel = 1
Gmodel = n_elements(GHI_arr)

d_clr = 35

CI = hici_hui(10^log_T_min)
RC = hiirca_hui(10^log_T_min)

xHI_noG = RC / (CI+RC)

xranges = dblarr(3,2)
yranges = dblarr(3,2)

yranges[0,*] = [-1.d-16,1.d-16]
xranges[0,*] = [-1.0d-2,1.0d-2]

yranges[1,*] = [-1.d-20,1.d-20]
xranges[1,*] = alog10([xHI_noG-0.5,xHI_noG+0.5])

yranges[2,*] = [-1.d-25,1.d-25]
xranges[2,*] = alog10([(1-1.0d-10),(1+1.0d-10)])


dxdt_no_tau0_arr = dblarr( 3, N, Tmodel, Gmodel )

xHI = dblarr( 3, N )

xHI[0,*] = dindgen(N)/(N-1) * (xranges[0,1] - xranges[0,0]) + xranges[0,0]
xHI[1,*] = dindgen(N)/(N-1) * (xranges[1,1] - xranges[1,0]) + xranges[1,0]
xHI[2,*] = dindgen(N)/(N-1) * (xranges[2,1] - xranges[2,0]) + xranges[2,0]
xHI = 10^xHI


if i_plot_multi then begin

    altay_set_x, windownum=0, xsize=900

    for i = 0, Tmodel-1 do begin
        T = 10^(log_T_min + i * d_log_T)
                
        for j = 0, Gmodel-1 do begin
            GHI = 10^(GHI_arr[j])
                     
            dxdt_no_tau0_arr[0,*,i,j] = dxdt_no_tau0( nH, T, GHI, y, xHI[0,*] )        
            dxdt_no_tau0_arr[1,*,i,j] = dxdt_no_tau0( nH, T, GHI, y, xHI[1,*] )        
            dxdt_no_tau0_arr[2,*,i,j] = dxdt_no_tau0( nH, T, GHI, y, xHI[2,*] ) 

        endfor
        
    endfor
    
    !P.multi=[0,3,2]

    ;---------------------
    altay_plot, [0],[0], yrange=yranges[0,*], $
      xrange=xranges[0,*], xtitle='log x_HI', ytitle='RHS', $
      charsize=3

    for j = 0, Gmodel-1 do begin
        altay_oplot, alog10(xHI[0,*]), dxdt_no_tau0_arr[0,*,0,j], thick=2, color=250-j*d_clr
    endfor
    
    oplot, [0,0], [-100,100], color=0
    oplot, [-100,100], [0,0], color=0

   
    ;---------------------
    altay_plot, [0],[0], yrange=yranges[1,*], $
      xrange=xranges[1,*], xtitle='log x_HI', ytitle='RHS', $
      charsize=3

    for j = 0, Gmodel-1 do begin
        altay_oplot, alog10(xHI[1,*]), dxdt_no_tau0_arr[1,*,0,j], $
          thick=2, color=250-j*d_clr
    endfor
    
    oplot, alog10( xHI_noG ) * [1,1], [-1,1], thick=2, linestyle=1, color=0

    oplot, [0,0], [-100,100], color=0
    oplot, [-100,100], [0,0], color=0

    ;---------------------
    altay_plot, [0],[0], yrange=yranges[2,*], $
      xrange=xranges[2,*], xtitle='log x_HI', ytitle='RHS', $
      charsize=3

    for j = 0, Gmodel-1 do begin
        altay_oplot, alog10(xHI[2,*]), dxdt_no_tau0_arr[2,*,0,j], $
          thick=2, color=250-j*d_clr
    endfor
    
    oplot, [0,0], [-100,100], color=0
    oplot, [-100,100], [0,0], color=0

endif





; analytic solution that considers particle contribution (xHI)
;==============================================================
;==============================================================
i_plot_multi = 1

x_floor = 1.0d-10
x_ceiln = 1.0d0 - x_floor

; for plotting several solutions
;------------------------------------------
x_tau_x_lo = dblarr(Tmodel,Gmodel)
x_tau_x_hi = dblarr(Tmodel,Gmodel)
tau0 = 1.0d1

dxdt_tau0_arr = dblarr( 3, N, Tmodel, Gmodel )

if i_plot_multi then begin


    for i = 0, Tmodel-1 do begin
        T = 10^(log_T_min + i * d_log_T)
                
        for j = 0, Gmodel-1 do begin
            GHI = 10^(GHI_arr[j])
            
            at_floor = dxdt_tau0_x_one( nH, T, GHI, tau0, y, x_floor )
            at_ceiln = dxdt_tau0_x_one( nH, T, GHI, tau0, y, x_ceiln )

            dxdt_tau0_arr[0,*,i,j] = dxdt_tau0_x_one( nH, T, GHI, tau0, y, xHI[0,*] )        
            dxdt_tau0_arr[1,*,i,j] = dxdt_tau0_x_one( nH, T, GHI, tau0, y, xHI[1,*] )        
            dxdt_tau0_arr[2,*,i,j] = dxdt_tau0_x_one( nH, T, GHI, tau0, y, xHI[2,*] ) 
                          
        endfor
        
    endfor

    ;---------------------
    altay_plot, [0],[0], yrange=yranges[0,*], $
      xrange=xranges[0,*], xtitle='log x_HI', ytitle='RHS', $
      charsize=3

    for j = 0, Gmodel-1 do begin
        altay_oplot, alog10(xHI[0,*]), dxdt_tau0_arr[0,*,0,j], thick=2, color=250-j*d_clr
    endfor
    
    oplot, [0,0], [-100,100], color=0
    oplot, [-100,100], [0,0], color=0

   
    ;---------------------
    altay_plot, [0],[0], yrange=yranges[1,*], $
      xrange=xranges[1,*], xtitle='log x_HI', ytitle='RHS', $
      charsize=3

    for j = 0, Gmodel-1 do begin
        altay_oplot, alog10(xHI[1,*]), dxdt_tau0_arr[1,*,0,j], $
          thick=2, color=250-j*d_clr
    endfor
    
    oplot, alog10( xHI_noG ) * [1,1], [-1,1], thick=2, linestyle=1, color=0

    oplot, [0,0], [-100,100], color=0
    oplot, [-100,100], [0,0], color=0

    ;---------------------
    altay_plot, [0],[0], yrange=yranges[2,*], $
      xrange=xranges[2,*], xtitle='log x_HI', ytitle='RHS', $
      charsize=3

    for j = 0, Gmodel-1 do begin
        altay_oplot, alog10(xHI[2,*]), dxdt_tau0_arr[2,*,0,j], $
          thick=2, color=250-j*d_clr
    endfor
    
    oplot, [0,0], [-100,100], color=0
    oplot, [-100,100], [0,0], color=0
    

endif

stop


GHI = 1.0d-12
T = 1.0d4
y = 0.0d0
nH = 1.0d-1

i_plot_multi = 1

x_floor = 1.0d-10
x_ceiln = 1.0d0 - x_floor

; for plotting several solutions
;------------------------------------------

N=100000
Tmodel = 5
Gmodel = 1
rhs = dblarr(N,Tmodel,Gmodel)

xHI = 10^( dindgen(N)/(N-1) * 11 - 9 )
tau0 = 1.0d1
i_plot_multi=1

if i_plot_multi then begin

    altay_set_x, windownum=2

    for i = 0, Tmodel-1 do begin
        T = 10^(4.0d0 + i * 0.1)
        
        print, 'T = ', T

        RC = hiirca_hui(T)
        CI = hici_hui(T)
        
        for j = 0, Gmodel-1 do begin
            GHI = 10^(-12 + j * 0.5)
            
            at_floor = ana_tau0_x_one( nH, T, GHI, tau0, y, x_floor )
            at_ceiln = ana_tau0_x_one( nH, T, GHI, tau0, y, x_ceiln )
            
            rhs[*,i,j] = ana_tau0_x_one( nH, T, GHI, tau0, y, xHI )        
            
        endfor
        
    endfor
    
    !P.multi=[0,1,3]

    ;---------------------
    altay_plot, [0],[0], yrange=[-1.d-13,3.d-14], $
      xrange=[-2.50,0.05]                        

    for i = 0, Tmodel-1 do begin
        altay_oplot, alog10(xHI), rhs[*,i,0], thick=2, color=250-i*40
    endfor
    
    oplot, [0,0], [-100,100], color=0
    oplot, [-100,100], [0,0], color=0

    ;---------------------
    altay_plot, [0],[0], yrange=[-1.d-15,1.d-15], $
      xrange=[-0.05,0.01]                        

    for i = 0, Tmodel-1 do begin
        altay_oplot, alog10(xHI), rhs[*,i,0], thick=2, color=250-i*40
    endfor
    
    oplot, [0,0], [-100,100], color=0
    oplot, [-100,100], [0,0], color=0

    ;---------------------
    altay_plot, [0],[0], yrange=[-1.d-16,1.d-16], $
      xrange=[-0.0001,0.0001]                        

    for i = 0, Tmodel-1 do begin
        altay_oplot, alog10(xHI), rhs[*,i,0], thick=2, color=250-i*40
    endfor
    
    oplot, [0,0], [-100,100], color=0
    oplot, [-100,100], [0,0], color=0

endif





; tau0 xHI << 1 - solving for xHI
;------------------------------------------------------------
R = GHI * tau0 + (RC + CI) * nH
Q = -( GHI + (CI+2*RC) * nH + (CI+RC)*nH*y )
P = RC * nH * (1+y)
d = Q*Q - 4 * P * R

xHIp = (-Q + sqrt(d) ) / (2*R)
xHIm = (-Q - sqrt(d) ) / (2*R)   ; this one

print, 'case of tau0 xHI << 1'
print
print, 'xHIp = ', xHIp, tau0 * xHIp
print, 'xHIm = ', xHIm, tau0 * xHIm
print

qq = - 0.5d0 * (Q + abs(Q)/Q * sqrt(d))

xHIp = qq / R
xHIm = P / qq   ; this one

print
print, 'xHIp = ', xHIp, tau0 * xHIp
print, 'xHIm = ', xHIm, tau0 * xHIm
print
print, '-----------------------------------------'

; NEGATIVE ROOT! 


; tau0 xHI >> 1 - solving for xHII
;------------------------------------------------------------
tau0 = 1.0d+2

R = -(GHI * tau0 * exp(-tau0) + (CI+RC) * nH)
Q = GHI * exp(-tau0) * (tau0-1) + CI * nH - (CI+RC) * nH * y 
P = GHI * exp(-tau0) + CI * nH * y
d = Q*Q - 4 * P * R


xHIp = (-Q + sqrt(d) ) / (2*R)
xHIm = (-Q - sqrt(d) ) / (2*R)   ; this one

print, 'caseof tau0 xHI >> 1'
print
print, 'xHIp = ', xHIp, tau0 * xHIp
print, 'xHIm = ', xHIm, tau0 * xHIm
print

qq = - 0.5d0 * (Q + abs(Q)/Q * sqrt(d))

xHIp = qq / R    ; this one
xHIm = P / qq

print
print, 'xHIp = ', xHIp, tau0 * xHIp
print, 'xHIm = ', xHIm, tau0 * xHIm
print
print, '-----------------------------------------'

; NEGATIVE ROOT! 





stop

; now try for single solution
;------------------------------------------


; set environment
;---------------------
nH = 1.0d-1
T = 1.0d4
GHI = 1.0d-12
tau0 = 1.0d1
y = 0.0d0



;zero_plot = ana_tau0_x_one( nH, T, GHI, tau0, y, xHI, dbg=1 )        
altay_plot, alog10(xHI), zero_plot, $
  yrange=[-1.d-14,1.d-14], /ys, $
  xrange=[-2,0.5], /xs, thick=1, color=250

altay_plot, alog10(xHI), zero_plot, thick=2, color=250, $
  yrange=[-1.d-13,1.d-13], /ys

altay_oplot, [-100,100], [0,0], color=0, thick=1
altay_oplot, [0,0], [-100,100], color=0, thick=1

stop

x_lo = 1.0d-8
x_hi = 1.0d0 - x_lo

stop

;at_zero = ana_tau0_x_one( nH, T, GHI, tau0, y, x_lo, dbg=1 )
;at_one  = ana_tau0_x_one( nH, T, GHI, tau0, y, x_hi, dbg=1 )
print, 'at zero = ', at_zero
print, 'at_one  = ', at_one

if (at_zero * at_one gt 0.0d0) then begin
    print, 'xHI = [0.0,1.0] bounds more or less than a single zero!'

    stop
endif


at_lo = at_zero
at_hi = at_one

; first check if we've bounded a zero by using the endpoints 
; xHI = x_lo and xHI = x_hi

if at_lo * at_hi ge 0.0d0 then begin
    print, ' bounded either 0 or more than one zero! '
    print, ' at lo: ', at_lo
    print, ' at hi: ', at_hi
endif


stop

; now we are sure to have the zero bounded

tol = 2.0d0
iclr=0
while tol gt 1.0d-4 do begin

    x_mid = 10^( (alog10(x_lo) + alog10(x_hi)) * 0.5d0 ) 
;    at_mid = ana_tau0_x_one( nH, T, GHI, tau0, y, x_mid, dbg=1 )
    
    
    altay_oplot, [alog10(x_mid),alog10(x_mid)], [-1,1], color=iclr

    print, 'lo/mid/hi: ', alog10(x_lo), alog10(x_mid), alog10(x_hi)
    print, 'lo/mid/hi: ', at_lo, at_mid, at_hi, tau0 * x_mid
    print

    if at_lo * at_mid gt 0 then begin
        x_lo  = x_mid
        at_lo = at_mid
    endif else begin
        x_hi  = x_mid
        at_hi = at_mid
    endelse
    
    tol = x_hi - x_lo

    iclr = iclr + 40

endwhile





end
